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Abstract 

A multi-cube method is developed for solving systems of elliptic and hyperbolic partial differ- 
ential equations numerically on manifolds with arbitrary spatial topologies. It is shown that any 
three-dimensional manifold can be represented as a set of non-overlapping cubic regions, plus a 
(~| ■ set of maps to identify the faces of adjoining regions. The differential structure on these mani- 

oc folds is fixed by specifying a smooth reference metric tensor Matching conditions that ensure the 

appropriate levels of continuity and differentiability across region boundaries are developed for 
arbitrary tensor fields. Standard numerical methods are then used to solve the equations with the 
appropriate boundary conditions, which are determined from these inter-region matching con- 
ditions. Numerical examples are presented which use pseudo-spectral methods to solve simple 
elliptic equations on multi-cube representations of manifolds with the topologies 
and S^. Examples are also presented of numerical solutions of simple hyperbolic equations on 
■ ^ . multi-cube manifolds with the topologies RxT^,RxS^ xS^ and RxS^. 

>>■ 
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1. Introduction 

The need to solve partial differential equations on manifolds having non-trivial spatial topolo- 
gies arises in many areas of physical science: from models of wormholes or the global structure 
of the universe in general relativity theory to global circulation models of the earth's atmo- 
sphere in meteorology and climatology. This paper develops practical methods for solving a 
variety of partial differential equations on manifolds having arbitrary spatial topologies. Every 
«-dimensional manifold (by definition) can be mapped locally into a portion of «-dimensional 
Euclidean space, R". A number of different numerical methods are capable of solving partial 
differential equations locally on open subsets of R". The topological structure of a manifold, 
however, affects the global solutions to partial differential equations in profound ways. This 
paper develops methods for fitting together local solutions, obtained from standard numerical 
methods, to form the desired global solutions on manifolds with arbitrary topologies. The dis- 
cussion here focuses on solving elliptic systems of equations on three-dimensional manifolds E 
with arbitrary topologies, and also hyperbolic systems of equations on four-dimensional mani- 
folds with topologies Rxl,. 

Solving partial differential equations numerically on manifolds with arbitrary topologies re- 
quires the creation of computational infrastructures (beyond those needed to solve the equations 
numerically on open subsets of R") that meet two basic requirements. The first requirement is 



Preprint submitted to Elsevier 



October 19, 2012 



that the manifold must be represented in a way that allows the points in the manifold, and the 
values of scalar and tensor fields defined at those points, to be referenced efficiently in a way 
that respects the underlying topological structure of the manifold. The second requirement is to 
create a way to specify the global differential structure of the manifold, i.e. the computational 
method must provide a way of representing globally continuous and differentiable scalar and 
tensor fields on these manifolds. The goal here is to develop practical methods that can be used 
on arbitrary manifolds by a wide range of different numerical methods. 

The first requirement is to find a systematic way of representing manifolds with arbitrary 
topologies. Every n-dimensional manifold can be mapped locally into a portion of n-dimensional 
Euclidean space R". For computational efficiency (and to avoid certain types of numerical in- 
stabilities) each manifold is represented here by a collection of non-overlapping n-dimensional 
cubes which cover the manifold, plus a set of maps that identify the faces of adjoining n-cubes. 
This decomposition is analogous to representing a manifold as a collection of non-intersecting 
«-simplexes (i.e., triangles for n - 2 and tetrahedrons for « - 3) that cover the manifold, plus 
maps that identify neighboring faces. Many numerical methods (including the pseudo-spectral 
methods used to produce illustrative examples for this paper) are easier to use in computational 
domains based on «-cubes rather than «-simplexes. Points in each of the n-cube regions are iden- 
tified by local Cartesian coordinates, and these coordinates are used to represent the solutions to 
the differential equations in each «-cube. This type of representation has been used for some time 
in numerical methods for solving partial differential equations on a two-sphere HH^jSIj and also 
in three-dimensional manifolds that are subsets ofB? 

Those ideas are gener- 
alized in Sec.|2] and it is shown that these generalizations can be applied to two-dimensional or 
three-dimensional manifolds having arbitrary topologies. Examples of these multi-cube repre- 
sentations are given in [Appendix A| for the three-dimensional manifolds with the topologies , 

The second requirement is to develop a method of representing (at least in the continuum 
limit) continuous and differentiable tensor fields on the multi-cube representations of manifolds 
developed in Sec. |2l Representing tensor fields within each of the «-cube regions is straightfor- 
ward: their components can be expressed in the tensor bases associated with the local Cartesian 
coordinates. These tensor components are functions of those coordinates, and their continuity (or 
differentiability) determines the continuity (or differentiability) of the tensor field itself. In gen- 
eral, however, the coordinate tensor bases associated with different «-cube regions are not even 
continuous (and can not be made continuous globally) across the interfaces that join them. The 
problem of defining the continuity and differentiability of tensor fields across n-cube interfaces 
is therefore non-trivial. The method introduced here makes use of a smooth reference metric 
tensor This reference metric must be supplied (along with the collection of n-cube regions and 
the associated interface maps) as part of the specification of a particular manifold. This metric 
is used to construct geometrical normal vectors at each interface, and these normals are used to 
construct the Jacobian matrices that map vectors (and tensors) across interfaces. The differentia- 
bility of tensors across the n-cube interfaces is defined in terms of the continuity of the covariant 
derivatives of those tensors, using the covariant derivative associated with the reference metric. 
The details of these continuity and differentiability conditions are given in Sec. [3] Examples 
of reference metrics which can be used to implement these continuity and differentiability con- 
ditions are given in Appendix A for the three-dimensional manifolds with the topologies , 

Systems of differential equations can be solved numerically on multi-cube representations 
of manifolds by fitting together local solutions from each n-cube region. The appropriate local 
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solutions are determined in each region by applying the correct boundary conditions on the n- 
cube faces. The appropriate boundary conditions are the ones that enforce the needed level of 
continuity and differentiability of the global solution at the region boundaries. These boundary 
conditions are developed in Sec.|4]for second-order strongly elliptic systems, and also for first- 
order symmetric hyperbolic systems of equations. These boundary conditions select the unique 
local solution in a particular n-cube that equals the desired global solution in that region. The 
collection of local solutions to the equations constructed in this way provides the desired global 
solution. 

The multi-cube method of solving systems of partial differential equations numerically on 
manifolds with non-trivial topologies is illustrated here by solving a series of test problems in 
Sees. |5] and |6] Simple second-order elliptic equations, and first-order symmetric hyperbolic 
equations, are solved numerically on manifolds with spatial topologies T^, x 5', and S^. 
These tests use pseudo-spectral methods to produce local solutions on each cubic region. The 
results are shown to converge exponentially (in an norm) to the exact global solutions (which 
are known analytically for these test problems) as the number of grid points used for the solution 
is increased. 



2. Building Multi-Cube Manifolds 

This section describes how «-dimensional manifolds can be represented using the multi-cube 
method. The idea is quite simple: n-dimensional multi-cube representations of manifolds consist 
of a set of non-overlapping n-cubes that cover the manifold, plus a set of maps that identify 
the boundary faces of neighboring cubes. An argument is presented in Sec. 12. II that all two- 
dimensional and all three-dimensional manifolds (with arbitrary topologies) can be represented 
in this way. A large class (but not all) higher-dimensional manifolds can also be represented using 
this multi-cube method. The multi-cube method provides a way of representing manifolds that 
facilitates the design of computational tools for solving partial differential equations on them. 
A simple infrastructure is introduced in Sec. 12.21 for systematically building, referencing and 
identifying the faces of the needed sets of «-cubes in these manifolds. These «-cube regions are 
joined together to form the desired topological manifold using maps that identify points on the 
faces of neighboring «-cubes. A simple framework for building and referencing these maps is 
presented. Only a small number of topologically distinct maps are needed for the case of three- 
dimensional manifolds (the main focus of this paper), and all of those maps are given explicitly. 



2.7. Existence of Multi-Cube Representations 

This subsection considers the question of whether two- and three-dimensional manifolds with 
arbitrary topologies admit multi-cube representations. The first step is to show that every two- 
manifold is homeomorphic to a set of squares (i.e. 2-cubes) glued together along their edges. 



The proof is based on the result of Rado 111 iL 11211 that all two-dimensional manifolds admit 



triangulations, i.e. that any two-manifold is homeomorphic to a set of triangles glued together 
along their edges. It is easy to show that a simple refinement of any triangulation on a two- 
dimensional manifold produces a multi-cube representation of that manifold. As illustrated in 
Fig lU let points "A", "B", and "C" denote the vertexes of one of the triangles in the triangulation. 
Add the midpoints of each edge of this triangle as additional vertexes, labeled "ab", "be", and 
"ac" in Fig. [T] Next, add the centroid of the triangle, the point labeled "d", and finally add 
as additional edges the line segments that connect "d" with the midpoints "ab", "be" and "ac". 
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The resulting complex consists of three quadrilaterals. When all of the triangles in a given 
triangulation are refined in this way, the result is a multi-cube representation of the two-manifold. 
The refinement consists of a set of quadrilaterals that are glued together edge to edge. Since the 
additional edge vertexes, "ab", etc. are always added at the geometrical midpoints, the edges 
of neighboring quadrilaterals constructed in this way will always coincide. These quadrilaterals 
are homeomorphic to squares (2-cubes). So the topological structure of a two-manifold can be 
thought of as a collection of non-overlapping 2-cubes that cover the manifold, plus a set of maps 
that identify the edges of adjoining 2-cubes. 



C . C 




Figure 1 : Each triangle in a triangulation of a two-dimensional manifold is refined by the addition of extra vertexes and 
edges to produce thi'ee quadrilaterals. This is done by first adding as new vertexes the midpoints of each edge, i.e. the 
points "ab", "be" and "ac" in the figure on the left. Next the centroid of the triangle, i.e. the point "d" in the figure on 
the right, is also added as a new vertex. Finally the line segments that join "d" to the midpoints "ab", "be", and "ac", the 
dashed lines rn the figure on the right, are added as new edges. 



A similar argument shows that every three-dimensional manifold has a multi-cube represen- 
tation, i.e. that every three-dimensional manifold is homeomorphic to a set of non-overlap ping 
"distorted" cubes glued together at their faces. The proof is based on a result of Moise lilAiOll 
that all three-dimensional manifolds admit triangulations by tetrahedrons, i.e. that any three- 
dimensional manifold is homeomorphic to a set of non-overlapping tetrahedrons glued together 
at their faces. It is easy to show that any tetrahedron can be decomposed into four "distorted" 
cubes glued together at their faces. (The term distorted cube is used here to describe a solid 
having six faces, each of which is a plane quadrilateral.) Distorted cubes are homeomorphic 
to geometrical cubes. It follows that every triangulation of a three-manifold can be refined (by 
adding appropriate vertexes, edges and faces) to obtain a multi-cube representation, i.e. a set of 
non-overlapping distorted cubes glued together at their faces. This argument demonstrates the 
existence of multi-cube representations for any three-dimensional manifold. 

The key to this argument is the representation of a single tetrahedron as four distorted cubes 
glued together. This can be done by refining the tetrahedron through the addition of vertexes, 
edges and faces as summarized in Fig. |2l Begin with a tetrahedron with vertexes labeled "A", 
"B", "C" and "D". First add vertexes to the midpoints of each edge, plus vertexes to the centroids 
of each face, the points "a", "b", "c" and "d" shown in the top left of Fig. |2l Adding the extra 
edges connecting "a", "b", "c" and "d" to the midpoints of each edge of the original tetrahedron 
completes the decomposition of each face into a set of distorted squares. Add one last vertex at 
the centroid of the tetrahedron, labeled "O" in the top right of Fig. |2] Connect "O" to the facial 
centroids, "a", "b", "c" and "d", by adding the edges shown as dash-dot line segments in the top 
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right of Fig. |2] Finally add the six internal quadrilateral faces that include the point "O" as an 
edge vertex. These additional vertexes, edges, and faces divide the tetrahedron into four volume 
regions (one adjacent to each tetrahedron vertex). The bottom of Fig.|2]shows these four regions 
more clearly. The regions adjacent to the vertexes "A" and "C" are shown with opaque faces, 
while those adjacent to "B" and "D" are shown with transparent faces. 




Figure 2: Top Left: Label the vertexes of the tetrahedron "A", "B", "C" and "D". Add vertexes at the midpoints of each 
edge, and additional vertexes at the centroid of each face of the tetrahedron, labeled "a" for the centroid of face "BCD", 
"b" for face "ACD", etc. Also add additional edges (shown as dashed line segments) connecting each centroid to the 
midpoint of each adjoining edge. Top Right: Add one additional vertex, labeled "O" at the centroid of the tetrahedron. 
Add additional edges (shown as dash-dot line segments) that connect "O" to the centroids of each face, and six additional 
faces that include "O" as a vertex. Bottom: Four "distorted" cubes that make up the tetrahedron are illustrated. The two 
cubes adjacent to vertexes "A" and "C" are shown with opaque shaded faces, while the faces of the cubes adjacent to "B" 
and "D" are transparent. 



Each of the four volume regions constructed above has six faces, and each of these faces has 
four edges and four vertexes. These faces are therefore quadrilaterals. It only remains to show 
that these quadrilaterals are planar. Call two edges of the original tetrahedron "complimentary" 
if they do not intersect at a vertex, e.g. the edges "AC" and "BD" are complimentary. Now 
consider the six bisecting planes of the tetrahedron, each one formed by an edge and the midpoint 
of the complementary edge of the tetrahedron. Each bisecting plane passes through the midpoint 
of the complementary edge, the centroid "O", as well as the facial centroids of the two faces 
adjacent to the complementary edge. For example, the bisecting plane formed by the edge "AC" 
and midpoint "bd" intersects "O" as well as the facial centroids "a" and "c". The quadrilateral 
formed by the vertexes "bd", "a", "O", and "c" is therefore a planar quadrilateral. It follows that 
each of the faces of the four volume regions is a planar quadrilateral, and therefore each volume 
region is a distorted cube. 

The vertexes added in this construction were placed at the geometric centroids of the trian- 
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gular faces, and at the centroid of the original tetrahedron. The edges added in this construction 
were also placed in geometrically determined ways: all of them along one of the bisecting planes 
of each edge of the original tetrahedron. These geometrically constructed features will therefore 
match on the triangular boundaries between neighboring tetrahedrons in any triangulation of a 
three-dimensional manifold. It follows that the distorted cubes constructed in this way will match 
face-to-face across all the tetrahedron boundaries as required for a multi-cube representation of 
the manifold. 

2.2. Infrastructure for Multi-Cube Manifolds 

Now turn to the problem of finding a systematic way of constructing multi-cube manifolds. 
The goal is to develop methods that can be used as part of the computational infrastructure 
for solving systems of partial differential equations on such manifolds. The discussion here 
is focused on three-dimensional manifolds S, but generalizations to other dimensions should 
be fairly straightforward. Let Sa denote a collection of geometrical cubic regions in R^. The 
subscript /I = {1, ...,A^) is used to label the individual regionsQ These cubes are used here as the 
domains of coordinate charts for the multi-cube representation of S. Let *Pa denote the invertible 
coordinate map that takes the region Sa into a subset of E: ^'a{Sa) c S. It will be useful to 
denote the boundary faces of these regions in as da^A, where a - +x denotes the faces 
intersecting the ±x axes, a - +y the faces intersecting the +y axes, etc. 

The discussion above shows that every three-manifold can be covered by a collection of 
non-overlapping cubes: \Ja^a{Sa) = 2. Non-overlapping here means that the images of the 
regions are non-intersecting, ^'a{Sa) n ^'b{Sb) - 0, for points in the interiors of Sa and Sb 
when A^B. It is convenient to choose the regions Sa in R^ to be scaled so they all have the 
same size L, and are all oriented along the same global Cartesian coordinate axes in R^. In this 
case the region Sa is completely determined therefore simply by specifying the location of its 
center ca = {c''^a,c^a,c"a) in R^- It is also convenient to arrange the regions Sa in R^ so they 
intersect (if at all) in R^ only at points on faces whose images also intersect in 2. In the multi- 
cube representations of manifolds satisfying these conditions, each point in the interior of the 
regions represents a unique point in 2, and each point in Z is the image of at least one point in the 
closure of UaSa- The Cartesian coordinates of R^ therefore provide a global way of identifying 
points in S. Tensor fields are represented on these multi-cube manifolds by giving the values of 
their components (expressed in the coordinate basis of R^) as functions of these global Cartesian 
coordinates. 

A multi-cube manifold consists of a set of cubic regions, Sa for a - {1, A^) that can be 
specified simply by giving the locations of their centers ca, along with a set of rules that deter- 
mine how the faces of these cubes are to be identified with one another When points on the 
images of two boundary faces *I'a(5„Sa) and ^BidpSB) intersect in Z, then the associated coor- 
dinate charts provide an invertible map from one boundary face to the other: daSA - ^^^idpSs) 
where 'I'^^ = ^F^' o ^b for points on the daSA and dpSB faces. Since the cubes Sa have uniform 
size and orientation in R^, there are only a small number of simple maps *P^^ needed to represent 
all the topologically distinct ways of mapping one face onto another. It is sufficient to consider 



' The term region in this paper is used to refer to the cubes Ba that form the basic topological structure of the manifold. 
It might be useful for computational efficiency to subdivide some (or all) of the cubic regions into a collection of smaller 
cubes, e.g. by cutting a cubic region into two, four, or eight smaller cubes. Those smaller cubic subsets of the Sa are 
referred to as subregions. 
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maps that identify the faces of two cubic region first by rigidly translating so the centers of the 
faces da^A and d/^Ss coincide, and then rigidly rotating and/or reflecting to align the two faces 
in the desired way. Thus it is sufficient to consider the simple maps that take the Cartesian 
coordinates x'^ of points in dpSB to the Cartesian coordinates x'^ of the corresponding points in 
daSA in the following way. 



r' - r' + f + C^'"(y' ~ r' - f'\ 
■^A ~ '-A ^ J a ^ ^Bp i^-^B J p ' 



(1) 



The vector ca + fa is the location of the center of the da^A face, and C^^ is the combined 
rotation and reflection matrix needed to achieve the desired orientation. Examples of the use of 



these methods is given in Appendix A where explicit multi-cube representations are constructed 
for manifolds with the topologies T \ S - x 5 ' and S ^ . 

Multi-cube manifolds are specified by giving the list of cubic regions Sa needed to cover the 
manifold, the vectors ca that determine the locations of their centers in R?, and the maps 
that determine how the regions are glued together. These maps, defined in Eq. ([1]), depend on the 
vectors ca and fa, and the matrix C^"p, so these quantities must all be specified to determine each 
map. The vector fa is the position of the center of the a face relative to the center of the region. 
Since the cubic regions are chosen to have uniform sizes and orientations, fa has the same form 
in each cubic region: 



A, = iL(± 1,0,0), 
fly - ^L(0,+1,0), 
/I, - ^L(0,0,+1), 



(2) 



where L is the size of the cubes. Since all of the cubic regions are aligned, the class of pos- 
sible rotations and reflections needed for Cg^ is quite small. These can all be constructed by 
combining 90-degree rotations about the normal to one of the faces, R^,, with mirror reflections 
about some (possibly different) direction. My?. Table [1] gives explicit expressions for the ma- 
trices that describe these elementary rotations and reflections in three dimensions. The most 
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general transformation of one face onto another can be constructed by taking products of these 
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elementary transformations. The group of possible C^^ in three dimensions generated in this 
way is therefore the octahedral symmetry group, (9/,, which has 48 distinct elements 11411 . The 
orientation preserving subgroup generated by the rotations alone has 24 elements. Note that 
Rq. ■ R-a = = = I, where I is the identity matrix. Since the number of possible maps 
constructed in this way is so small, it is easy to write a flexible code that is capable of setting 
up the multi-cube structures and all the needed gluing maps for three-manifolds with arbitrary 
topologies. 



3. Specifying Differential Structures on Multi-Cube Manifolds 

This section describes a practical and efficient way to define C'' differential structures on 
multi-cube manifolds. It is useful to begin with a brief discussion of the traditional way such 
structures are defined. The differential structure on a manifold provides the framework needed 
to represent differentiable scalar and tensor fields on that manifold. The usual method of spec- 
ifying a differential structure is to cover the manifold with a set of overlapping domains Da, 
and set of maps that assign coordinates to the points in each domain: T^'(2)^) c R". These 
coordinate maps provide a differential structure for the manifold if they have the property that 
the composition maps Tg = T^' o are differentiable (or C*^^') transformations from the co- 
ordinates of one patch to the other for points in the overlap Da n Db- The Jacobian matrices 
associated with these coordinate transformations J'^'j = dx'^/dxg determine the transformations 
for C* differentiable tensors from one coordinate representation to another in these overlaps. 

It is possible to use the traditional method of defining differential structures on multi-cube 
manifolds, but to do so requires that non-trivial additional structures must be added to the basic 
multi-cube construction (since the domains that define that basic structure do not overlap). The 
most straightforward approach would be to require that each multi-cube manifold be provided 
with an additional set of overlapping domains Da ^ ^a{Sa) and a set of C*^' related coordinate 
maps for the new overlapping Da domains. An alternative, more minimalist, approach would 
be to require that suitable Jacobian matrices 7^^", in addition to the connection maps 'P^^, be 
provided on each interface between regions in multi-cube manifolds. This minimal structure 
would provide the transformations needed to define differentiable scalar and continuous tensor 
fields on these manifolds. If C*^' differentiable scalars or differentiable tensor fields are 
needed, then in addition to 7^^", all of their ^* order derivatives dgJ'^'j would also have to be 
specified on each interface between regions. 

It might seem redundant and unnecessary to require that the Jacobian matrices J^^'j and their 
derivatives be specified on the interfaces in multi-cube manifolds, in addition to the interface 
coordinate maps *I'^j^ defined in Eq. ([1]). After all, the Jacobian matrices associated with those 
interface maps, J^'^'j - C^^'j, and their derivatives, dskJefj - dskCg^'j = 0, could be used to trans- 
form tensor fields at the boundary interfaces. Unfortunately it is easy to see that the coordinate 
maps ^'a used in Sec. |2] to construct the multi-cubes are not suitable for constructing a global 
C'' differential structure on most manifolds. If they were, the basis vectors Sai associated with 
these coordinates would be smooth global non-vanishing vector fields. These vector fields could 
be used in this case to construct a global smooth flat metric on the manifold. Since most mani- 
folds do not admit global flat metrics, the existence of a complete set of smooth non-vanishing 
coordinate vector fields can not exist on most manifolds. Figure [3] drawn from the perspective 
of a smooth coordinate patch that covers both sides of an interface boundary, illustrates how 
the multi-cube coordinates in neighboring regions can be continuous while failing to be differ- 
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entiable across region boundaries. The coordinate region Si on the left, matches to coordinate 
region S2 on the right across the Xi = X2 interface in Fig. [3] The coordinate vectors tangent to 
this interface, e.g. dr^ and Sy^, are continuous across this interface, while those not tangent to the 
boundary, i.e. 5x, and are discontinuous there. 




Figure 3: Maps 'Va define continuous but (typically) non-differentiable transitions between cubic regions. This example 
shows that the basis vectors tangent to the boundary, dy^ and 5)', , are continuous, while those not tangent to the boundary, 
and > ai'e not. 

Both approaches described above for specifying differential structures on a multi-cube man- 
ifolds require that a great deal of extra structure be provided. This paper proposes a third, more 
elegant and more efficient, approach that can be incorporated more easily into the computational 
infrastructure for solving partial differential equations numerically. Every manifold with a C*^^' 
differential structure admits a symmetric positive definite C*^ differentiable metric tensor gij. The 
method proposed here for specifying the global differential structure on a multi-cube manifold 
requires that the components of (any) one of these diff'erentiable reference metrics, gij, be 
provided in the global Cartesian coordinate basis used to define the multi-cube manifold. The 
components of this reference metric gij will be C* functions of the multi-cube Cartesian coordi- 
nates within each region Sa, but will (in general) be discontinuous across the interfaces between 
regions. The only requirement on this reference metric is that it must be sufficiently smooth, C*, 
when represented in a global C*^' coordinate atlas. The C*^' coordinate charts themselves 
need not be given as part of the specification of the multi-cube manifold. Their only use in this 
method is to ensure a priori that the reference metric meets the needed smoothness requirements. 

Once a suitable reference metric gij is provided, it is straightforward to construct the Jaco- 
bian matrices Ji'^'- and the dual Jacobian matrices J*f^^ needed to transform continuous tensor 

Dp J Aai 

fields across the interface boundaries in multi-cube manifolds. Assume that the da^A bound- 
ary of region Sa is identified with the d/^SB boundary of region Sb by the map ^F^^ given in 
Eq. ([1]). The transformation taking the region Sb representation of a vector into the region Sa 
representation at one of these identified boundary points is an expression of the form 



where 7^"'. is in effect the Jacobian matrix of the transformation. The analogous transformation 
law for covectors wbi is. 



I _ jAai j 



(3) 




(4) 



where J' 



Aai 



is in efi'ect the dual Jacobian matrix. 
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Let gAij denote the coordinate components of the reference metric in the muhi-cube coordi- 
nate basis of region Sa, and let tiAoi denote the outward directed normal covector to the surface 
daSA- This interface is a surface of constant coordinate jcj^, so the geometrical normal covector 
is proportional to (9/i,x^. The normal covector is therefore given by 

«Am - I (5) 
^]8A9Aj^^AkX'^ 

where g'^ is the inverse of the reference metric gAij- The sign is chosen in this expression to make 
iiAai the outgoing unit normal. The unit normal vector n'^^ is related to riAai by n'^^ = g'^yiAaj.- 

The Jacobian matrices needed to transform vectors and covectors (and therefore any type 
of tensor field) across boundary interfaces are simple functions of the quantities C^^'j and C^J. 
(which define the identification maps ^^"p), as well as the normals to the boundary surface, n^^, 

■1%) = Cg"p[[5)-n)^pnBii^-n\^nBpj, (6) 

Jfai' = ((^f - Cff^ - "A^riW^^. (7) 

The Jacobian matrices defined in Eqs. ^ and (|7]i are the unique ones with the properties: 
a) They map the geometrical normals into and nspj into -iiAai, 

n'Aa = -JbpHp, (8) 

riAai = -JlT^BBj, (9) 



(i.e. the outward directed normal of one region is identified with the inward directed normal of 
its neighbor). 

b) The Jacobian matrix J^^'j transforms any vector t' tangent to the boundary (i.e. any vector 

' 8/3 



satisfying f'n, - 0) using the continuity of the ^PoJ? maps: 



c) The Jacobian matrix Jg^j and its dual J*^^' are inverses 

cAi jAtyi T*Bpk / 1 1 \ 

^Aj^ JBpkhaj ■ (11) 

This last property ensures that tensor contractions and traces transform properly under these 
boundary interface mappings. 

The Jacobian matrices constructed in Eqs. Q and (|7]) using the identification maps '^'^^ and 
the reference metric gij define the transformations needed to connect arbitrary tensor fields across 
the interface boundaries of multi-cube manifolds. These transformations make it possible there- 
fore to define what it means for a global tensor field to be continuous on multi-cube manifolds: A 
tensor field is continuous on a multi-cube manifold if its multi-cube coordinate components are 
continuous within each region 3a, and if its multi-cube coordinate components at each interface 
boundary point are equal to the transform of its components from the neighboring region. 
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The reference metric can also be used to define a smooth connection 



(12) 



that can be used to define a covariant derivative operator V,-. This covariant derivative is related 
to the coordinate partial derivatives (within each region Sa) by the usual expressions for the case 
of vectors and covectors: 



The covariant gradients of tensors, e.g. V,v-' and V,Wj, are themselves tensor fields. Therefore 
they are transformed at interface boundaries using the Jacobian matrices defined in Eqs. (|6]l and 
(|7]i as well. Thus, for example, the gradients of vectors and covectors transform as. 



Using these transformation laws it is straightforward to define what it means for a global tensor 
field to be differentiable on a multi-cube manifold: A tensor field is differentiable if the tensor 
and its covariant gradient are continuous everywhere including across all multi-cube interfaces. 
The concept of C* tensors can be built up in a straightforward way simply by taking A:* order 
covariant gradients of tensors and demanding that the tensor and all gradients up through A:* 
order be continuous global tensor fields. 

The addition of a smooth (i.e. C*^ differentiable) positive definite reference metric gij there- 
fore provides all the additional information needed to define a global differential structure on 
any multi-cube manifold. 

4. Interface Boundary Conditions for Multi-Cube Manifolds 

The multi-cube representations of manifolds provide a practical framework in which to solve 
systems of partial differential equations numerically on manifolds with non-trivial spatial topolo- 
gies. The idea is to solve those equations on each of the cubic regions Sa separately, using 
boundary conditions on the faces that ensure the combination of local solutions from each 
region satisfies the system of equations globally — including at the boundaries. Solving differ- 
ential equations using multi-patch methods is a common practice in computational physics on 
manifolds that are subsets of Such methods are used for example in 

the pseudo-spectral code SpEC (developed by the Caltech/Cornell numerical relativity collabo- 
ration [15, 16, 17, 18, il^) to solve Einstein's equations. The multi-cube framework developed 
here extends the class of problems accessible to such codes by allowing them to solve prob- 
lems on computational domains that can not be covered by a single global coordinate chart. 
This generalization provides a method of solving differential equations on two-dimensional and 
three-dimensional manifolds with arbitrary topologies, in addition to a very large class of higher 
dimensional manifolds. The code changes needed to implement these more general multi-cube 
methods require fairly minor generalizations of the way boundary conditions are imposed at 
the interfaces between cubic regions in standard multi-patch codes. The needed generalizations 
are described here in some detail for second-order quasi-linear strongly-elliptic and first-order 
symmetric-hyperbolic systems of equations. 



(13) 
(14) 




(15) 



(16) 
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4.1. Interface Boundary Conditions for Elliptic Systems 

A second-order quasi-linear strongly-elliptic system of equations for a collection of tensor 
fields ic^ can be written in the form 

V; [M-'*-^s(u)VtM®] = F'^ivL, Vu), (17) 

where V, is some covariant derivative operator, M^*^s(u) may depend on the fields but not their 
derivatives, and F®(u, Vu) may depend on the fields and their first derivatives. The script indexes 
'S,c, ... in these expressions label the components of the collection of tensor fields that make 
up u^. Such a system is strongly elliptic if there is a positive definite metric on the space of 
fields, Sj['B, a positive definite spatial metric, g'^, on the manifold (e.g. the reference metric used 
to define the multi-cube structure) and a positive constant, C > 0, such that 

w,w,M-'%5c» v^v® > Cg^'^WjWk Syts v^^P (18) 

for every and every Wj lEoll . 

All differentiable soltuions to second-order elliptic systems of this type are smooth, assum- 
ing the quantities M-'^^^g and are smooth ll20ll . Boundary conditions for these equations at 
internal inter-region boundaries are therefore quite simple: the solutions ir^ and their normal 
derivatives n'Viir^ (where «' is the normal to the boundary) must be continuous when trans- 
formed appropriately across inter-region boundaries. 

These continuity conditions can only be imposed at the interface boundaries by transforming 
the fields computed in one region, Sb, into the tensor basis used by its neighboring region, Sa . 
The fields are (by assumption) a collection of tensor fields whose components are transformed 
across region boundaries using the Jacobian as defined in Eqs. (O and (|4]i. Thus the fields 
(expressed in the tensor basis associated with the coordinates x'g from the region Sb) are related 
to the fields (in the tensor basis associated with the coordinates x'^ from the region Sa) by a 
transformation of the form, 

^^J-^swf, (19) 

where J''% is the multi-component Jacobian appropriate for each tensor part of u^. For exam- 
ple, a system whose fields consist of a scalar, a vector, and a covector m® = {t//,v',Wi}, would 
transform as follows, 

J^s4-[>f^s, J'^V^B^ CwB,}. (20) 

The boundary conditions for second-order elliptic systems also place conditions on the nor- 
mal derivatives of the fields, n'ViU^. The covariant gradient of a tensor field is itself a tensor 
field, so these gradients are transformed across region boundaries by an equation analogous to 
Eq. (HI: 

Va/"a - JAai J'B ^BjUb- (/i) 

It may be more convenient in some cases to impose the needed continuity conditions on the par- 
tial derivatives, n'diU^, rather than the covariant derivatives of the fields, n'Vtu^. The interface 
boundary transformations needed in this case are easy to obtain from Eq. ( ISTT i: the covariant 
derivatives ^Ak and V^^: that appear in this condition are re-expressed in terms of the partial 
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derivatives Oa; and Bm, and the connection coefficients F^^.^ and T^^.^. For the case of vector and 
co-vector fields, the resuhing partial derivative transformation laws are given by, 

OAkVA - JAak hpi OB(Vg + [J ^^j^ J b/},, 1 BC j " J B/i j ^ Aknj^B' (^^) 
OAkWAi - J^ak •'Aai ^BCVBj ' [J a^,, J Ad ^ Bin ' J Aan ^ Aki) ^Bj- 

The needed interface boundary conditions for second-order elliptic systems can now be stated 
precisely: Let Sa and Sb represent cubic regions whose faces da^A and d/^SB are identified. Let 

and denote the fields evaluated in the cubic regions Sa and Sb respectively. The required 
interface boundary conditions can then be written as, 

«f=J''5i<, (24) 
to be imposed on the boundary face dpSs, and the equation, 

n\^AiU^ = n'^rlfj^'B^Bkul, (25) 

to be imposed on the boundary face da^A- 

The required continuity conditions can be imposed numerically by replacing the elliptic sys- 
tem, Eq. (T7\ . with the equation for the continuity of the fields on the grid points of one of the 
boundary faces, S/^Sb, and the equation for the continuity of the normal derivatives on the grid 
points of the other face Together these boundary conditions ensure that the global solution 

to Eq. ( fTTl i will have the required continuity and differentiability at interface boundaries. Second- 
order strongly-elliptic systems can be solved using either Dirichlet or Neumann type boundary 
conditions. Thus the continuity conditions imposed here are exactly those needed to ensure the 
well-posedness of the boundary value problem within each cubic region. 

Boundary conditions of this type are already used successfully and routinely in elliptic-solver 
codes that implement traditional multi-patch methods (see e.g. Ref. |fl6[l). The only difference 
between the boundary conditions used in those traditional multi-patch codes and the ones in- 
troduced here is the form of the Jacobian matrices used to transform the components of tensors 
and their derivatives at the interfaces between regions. In traditional multi-patch methods these 
Jacobians are just identity matrices, because in those cases there was always a smooth global 
coordinate basis that could be used to represent tensor fields in all computational subdomains. 
In the multi-cube method introduced here, these Jacobians contain critical information about the 
differential topology of the manifold. 



4.2. Interface Boundary Conditions for Hyperbolic Systems 

A first-order symmetric-hyperbolic system of equations for the dynamical fields (assumed 
here to be a collection of tensor fields) can be written in the form 

<9,M-^ + A*-^s(u) V^M® = F'^iu), (26) 

where the characteristic matrix, A'^-^s(u), and source term, F^(u), may depend on the fields ir^ 
but not their derivatives. The script indexes j{, s, c, ... in these expressions label the components 
of the collection of tensor fields that make up ic^. These systems are called symmetric because, 
by assumption, there exists a positive definite metric on the space of fields, Sjis, that can be used 
to transform the characteristic matrix into a symmetric form: SjicA'^'^^ = A* - A* . 
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Boundary conditions for symmetric -hyperbolic systems must be imposed on the incoming 
characteristic fields of the system. The characteristic fields M^(whose index ■k labels the collec- 
tion of characteristic fields) are projections of the dynamical fields onto the left eigenvectors 



of the characteristic matrix (cf. Refs. ||21U22|1 ). 
defined by the equation. 



= e'^y^('a) m^, (27) 



e%(n) rikA'^^siu) = e%(n). (28) 

The co-vector «i that appears in this definition is the outward pointing unit normal to the surface 
on which the characteristic fields are evaluated. The eigenvalues \'cx) are often referred to as the 
characteristic speeds of the system. The characteristic fields represent the independent dy- 
namical degrees of freedom at the boundaries. These characteristic fields propagate at the speeds 
V{-K) (in the short wavelength limit), so boundary conditions must be given for each incoming 
characteristic field, i.e., for each field with speed V(;k) < 0. No boundary condition is required (or 
allowed) for outgoing characteristic fields, i.e., for any field with v^c, > 0. 

The boundary conditions on the dynamical fields that ensure the equations are satisfied 
across the faces of adjoining cubic regions are quite simple: data for the incoming characteristic 
fields at the boundary of one region are supplied by the outgoing characteristic fields from the 
neighboring region. The boundary conditions at an interface between cubic regions require that 
the dynamical fields in region !Ba be transformed into the tensor basis used in the neighboring 
region Sb- When the dynamical fields are a collection of tensor fields (as assumed here) their 
components are transformed from one coordinate representation to another using the Jacobian of 
the transformation as described in Eq. ( fT9l l. In this case the needed boundary conditions can be 
stated precisely for hyperbolic evolution problems: Consider two cubic regions Sa and Sb whose 
boundaries and dp^B are identified by the map ^F^^ as defined in Eq. ([1]). The required 

boundary conditions on the dynamical fields consist of fixing the incoming characteristic 
fields mJ, i.e., those with speeds v^k) < 0, at the boundary da^A with data, m^, from the fields on 
the neighboring boundary S/^Sb- 

u'^A - e'^^irDJ^ul (29) 

The matrix of eigenvectors, e'^yiin), that appears in this expression is to be evaluated using the 
fields from region Sb that have been transformed into region Sa where the boundary condition is 
to be imposed. This boundary condition must be applied to each incoming characteristic field on 
each internal cube face, i.e., on each face that is identified with the face of a neighboring region. 

This type of boundary condition is used routinely and successfully by hyperbolic evolution 
codes, such as the Caltech/Cornell SpEC code, that implement traditional multi-patch methods. 
Those traditional apphcations differ from the multi-cube methods discussed here only in the fact 
that tensors in those traditional cases could always be expressed in terms of the global coordinate 
basis. The generalized Jacobians JT^s needed to transform tensors across interface boundaries 
in those traditional applications of multi-patch methods are therefore just the identity map. In the 
more general multi-cube construction introduced in Secs.|2]and[3l the Jacobians contain critical 
information about the differential topology of the manifold, so the transformations used here 
must be slightly more complicated than those used in the traditional multi-patch case. Other than 
that simple difference, however, the boundary conditions introduced here are the same as those 
used in the traditional multi-patch methods. 
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In some cases, like systems representing second-order tensor wave equations, the dynami- 
cal fields will include a collection of primary tensor fields plus a collection of secondary fields 
representing the first derivatives of the primary fields. In most cases the secondary fields can 
be defined using a covariant derivative, thus making them tensor fields as well. The Einstein 
equations are somewhat problematic, because the most natural covariant derivative of the metric 
tensor (the primary tensor field in this case) vanishes identically. Thus first-order symmetric- 



hyperbolic representations of the Einstein equations are not generally co-variant 112211 . They can 
be made fully covariant however by defining the secondary dynamical fields using the covariant 
derivative associated with the non-dynamical reference metric that defines the differential topol- 
ogy of the manifold. This type of fully covariant first-order representation of the Einstein system 
will be discussed in detail in a future publication. 

5. Numerical Tests of a Multi-Cube Elliptic Equation Solver 

This section discusses a series of tests of the numerical solution of elliptic equations on 
compact three-manifolds using the multi-cube methods described in Sees. |2] [3] andH) These 
tests find numerical solutions to the equation 

V'V,iA - cV = /, (30) 

where tf) is a scalar field, V, represents the covariant derivative associated with a fixed smooth 
positive-definite metric gij on a particular three-manifold, c is a constant, and / is a fixed source 
function. The constant term, with c- > 0, ensures the solution to this equation is unique on 
any compact three-manifold. This equation is solved here on the three-manifolds whose multi- 
cube representations are described in Appendix A with a flat metric, 5 ^ x 5 ' with a round 



constant-curvature metric, and S ^ with the standard round constant-curvature metric. The source 
functions / for these tests are chosen to ensure that the solutions are non-trivial functions 
which are known analytically. 

The accuracy and effectiveness of the numerical solutions of Eq. ( l30t are evaluated in two 
ways. The first accuracy indicator used here is the residual, R, which measures how well the 
numerical solutions satisfy the discrete form of the differential equations. This numerical residual 
is defined as 

= V'V,iAiV-cViv-/, (31) 

where i^ai is the numerical solution of the discrete form of Eq. ( l30l l. The size of this residual is 
monitored for each numerical solution by evaluating its norm and computing the normalized 
residual error quantity, £«, defined as 



A -r ■ (32) 



The second accuracy indicator used here measures the error in the numerical solution itself: 
Aif/ - ^E—^N, where ijjE and represent the exact analytical solution and the discrete numerical 
solutions respectively. The magnitude of Ai/^ is evaluated using the scale invariant measure of 
the solution error: 



/(A^)^Vg^3^ 

A r ^ ■ \->->) 
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The numerical tests described here were performed using the elliptic equation solver that is 
part of the SpEC code lfl6ll . This code, developed originally by the Caltech/Cornell numerical 
relativity collaboration, uses pseudo-spectral methods to represent functions and evaluate their 
spatial derivatives. It solves elliptic equations using the PETSc toolkit of linear and non-linear 
equation solvers. Each cubic region in the tests described here is subdivided into one or more 
computational subregions, on which field components are represented using Chebyshev basis 
functions at the Gauss-Lobatto collocation points. The structure of these subregions was chosen 
to achieve fairly uniform spatial resolution. The particular choice of subregions is described in 
the discussion of each test. 

These numerical tests verify that several new ideas introduced in Secs.[2l[3l|4land [Appendix A 



are correct, and that these ideas have been implemented correctly in the SpEC code. The most 
fundamental new ideas tested here are the inter-region boundary conditions, Eqs. (l24l) and dZSl ), 
for elliptic equations. These internal boundary conditions depend on the Jacobians and their 
derivatives, which depend in turn on the inter-region boundary maps in a critical way for man- 
ifolds with non-trivial topologies. These Jacobian terms contribute to the boundary conditions 
in a non-trivial way even for the simple scalar elliptic equation ( l30b used in these tests. These 
tests also depend in a non-trivial way on the multi-cube representations of the reference metrics 
Eqs. ( IA.9I ) and ( IA.20l i and their associated covariant derivatives on the manifolds 5^ x 5 ' and 
5^. If any of these new elements of the multi-cube method were incorrect (or were implemented 
incorrectly in the code) the numerical tests described here would not achieve the exponential 
convergence in the solution error measure &^ that is seen in these tests. 

5.7. Tests of a Multi-Cube Elliptic Equation Solver on 

The numerical tests described here use the multi-cube representation of the three-manifold 
with topology given in Appendix A. 1 The reference metric in this case is the flat Euclidean 



metric, Eq. ( lA.ll ). so the covariant derivatives which appear in the elliptic Eq. ( l30l l are just the 
Cartesian coordinate partial derivatives. When written in terms of the multi-cube Cartesian co- 
ordinates on T^, therefore, this equation takes the simple form, 

V'V,iA - cV = dl(f/ + d^iff + 5^i/r - cV = /■ (34) 

This equation is solved numerically in these tests using the source function / given by. 



fix,y,z) = -{cj^ + c^)cos 



In 

— (kx + £y + mz) 



(35) 



where k, £, and m are integers, c is a constant c - 1/L, and co is given by 

\2 



J^j (;tHf2^m2). (36) 



The exact analytical solution to this equation is given by 

'2jt 



^Eix,y,z) = cos 



^ (kx + £y + mz) 



(37) 



The numerical tests of the solutions to Eqs. (I34li-(l36ll were performed using a source function 
with k - { - m — 2. These tests were performed on a set of eight computational subregions 
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using a range of numerical resolutions having N - S, 10, 12, 14, 16, 18 and 20 collocation points 
respectively in each spatial direction in each subregion. These subregions divide the one cubic 
region Si needed to represent 7"-^ into eight cubes: each half the size of the region in each spatial 
direction. The internal boundary maps between these subregions are just the trivial identity 
maps. The graphs of the solution errors and the residual errors as defined in Eqs. (l32T l 
and ( l33T l, for these tests are shown in Fig. |4] The elliptic sover for these tests were run until the 
residual errors Sn were reduced to the level of numerical roundoff. These results demonstrate that 
the boundary conditions introduced here on region boundaries were implemented correctly and 
efficiently: successfully achieving the exponential convergence expected of spectral numerical 
methods. 




Figure 4: Errors in the numerical solutions Ai// of the elliptic Eq. )34t on with k 



m = 2, as quantified by 



the error measures and £«. The parameter N is the number of collocation points used for these tests in each spatial 
direction in each computational subregion. 



5.2. Tests of a Multi-Cube Elliptic Equation Solver onS^xS^ 

The numerical tests described here use the multi-cube representation of the three-manifold 
with topology S' x given in Appendix A. 2 The reference metric used in this case is the 



constant-curvature round metric given in terms of angular coordinates {x, 0, (p) in Eq. ( IA.8I ). and 
in the multi-cube Cartesian coordinates used in these tests in Eq. iA.9i . This choice of reference 
metric makes the elliptic Eq. ( l30t somewhat more complicated in this case. In terms of the 
standard angular coordinates this equation has the form 



2 , _ ^ dg [singgfllA] ^l^ _ ^2.. 



V'V,,A-cV- 4r + 7 + „, 9 -^V-/- (38) 



This equation is solved numerically in these tests with a source function / given by, 

fix, 0, if) = -{or + c'-)% [e''^YiJ6, ^)] , (39) 

where Ycm(0,tf) is the standard S~ spherical harmonic function, k, £, and m are integers, c is a 
constant c - I/R2, to is given by 

, e{e+i) k- 

OJ^ = , + ^, (40 
Rl R\ 
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and %[Q] denotes the real part of a quantity Q. The exact analytical solution to this equation is 
given by 



(41) 



The numerical solution to this equation is carried out using the Cartesian coordinates of the 
multi-cube description of 5^ x 5 ' described in Appendix A. 2 The covariant derivatives used 



by the SpEC code for this test are evaluated using the Cartesian coordinate representation of the 
round metric given in Eq. iA.9i . The source function / that appears on the right side of Eq. ([38]), 
is evaluated in the multi-cube Cartesian coordinates used for these tests with the transformations 
between the angular and Cartesian coordinates given in Tables lA~4l and lA3] 

The tests performed here used the source function given in Eqs. ([39])-(l40t with k - ( - m - 
2. These tests used a set of twelve computational subregions to represent the six cubic regions of 
5^ X 5 cf. Fig. IA.10l These subregions divide each region in the periodically identified z direc- 
tion into two subregions. These tests were performed using N 10, 12, 14, 16, 18, 20 and 22 
collocation points respectively in each spatial direction in each of the computational subregions. 
The boundary conditions at the inter-region boundaries are based on the maps specified in Ta- 
ble |A3] The graphs of the solution errors £^ and the residual errors as defined in Eqs. ( l32l i 
and ( |33] ). for these tests are shown in Fig.|5] The elliptic sover for these tests were run until the 
residual errors &n were reduced to the level of numerical roundoff. This graph demonstrates, for 
the non-trivial 5 ~ x 5 ' case, that the computational region boundary conditions developed here 
have been implemented correctly and efficiently, achieving the exponential convergence expected 
of spectral numerical methods. 
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Figure 5; Errors in the numerical solutions Ai/< of the elliptic Eq. )38t on 5^ X S ' with k = t = m = 2, quantified by 
the eiTor measures £^ and £r. The parameter N is the number of collocation points used for these tests in each spatial 
direction in each computational subregion. 



5.3. Tests of a Multi-Cube Elliptic Equation Solver on S 

The numerical tests described here use the multi-cube representation of the three-manifold 
with topology given in Appendix A. 3 The reference metric used in this case is the stan- 
dard constant-curvature round metric for S ^ given in terms of angular coordinates {x, 0, <p} in 
Eq. ( IA.19I ). and in the multi-cube Cartesian coordinates used in these tests in Eq. ( IA.20I I. This 
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choice of reference metric fixes the elliptic Eq. ( l30l ) to have the form, 



V'V,)/' - cV 



5^[sin2^5^^] 54sin05fliA] 



R\5\vl X RlsmOsivrx ^? sin sin ;t' 



(42) 



when expressed in terms of the standard angular coordinates {x,6,(f] used on 5^. The source 
function / used in these numerical tests is given by. 



(43) 



where the Yumix^ ^re the spherical harmonics described in [Appendix B k, i, and m are 
integers, c is a constant c - \ /Rt,, and a> is given by 



, k{k + 2) 

W - T . 

Rl 



The exact analytical solution to this equation is given by 



(44) 



(45) 



The numerical solutions of Eq. (l42l i are carried out for these tests using the multi-cube rep- 
resentation of described in Appendix A. 3 The covariant derivatives used by the SpEC code 



for this test are evaluated using the multi-cube Cartesian coordinate representation of the round 
metric on 5^ given in Eq. ( IA.20I ). The source function /, defined in Eq. ( |43] ), is evaluated in 
terms of the multi-cube Cartesian coordinates for these tests using the transformations between 
the angular and the Cartesian coordinates given in Tables lA!8] and |A.9l 




Figure 6: Errors in the numerical solutions At// of the elliptic Eq. i42i on S ^ with i = f = m = 2 as quantified by the error 
measures £^ and Sr. The parameter N is the number of collocation points used for these tests in each spatial direction 
in each computational subregion. 



The numerical tests described here solved the elliptic Eqs. (I42l) - (l44l) with the parameter val- 
ues k - £ = m = 2 in the source function /. These tests were done using a set of eight 
computational subregions, corresponding to the eight cubic regions needed to represent S^, cf. 
Fig. lA. Ill These tests used = 8, 10, 12, 14, 16, 18, 20 and 22 collocation points respectively 
in each spatial direction in each of the computational subregions. The boundary conditions at 
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the region boundaries for these tests are based on the interface identification maps specified in 
Table lA. 8] The graphs of the solution errors &^ and the residual errors &r, defined in Eqs. ( |32] | 
and ( 1331 ). for these tests are shown in Fig.|6] The elliptic sover for these tests were run until the 
residual errors &r were reduced to the level of numerical roundoff. This graph demonstrates for 
another non-trivial example that the inter-region boundary conditions developed here have been 
implemented correctly and efficiently. Figure|6]also demonstrates that these numerical tests have 
achieved the exponential convergence expected of spectral numerical methods. 

6. Numerical Tests of a Multi-Cube Hyperbolic Equation Solver 

This section discusses numerical tests of the multi-cube methods for solving hyperbolic evo- 
lution equations on compact three-manifolds as described in Secs.|2]|3] and|4] These tests find 
numerical solutions to the scalar wave equation 

- (9,V + V'V,iA = 0, (46) 

where V, represents the spatial covariant derivative on the fixed geometry of the spatial three- 
manifold. This equation is solved here on the three-manifolds described in [Appendix A[ 
with a flat metric, 5 ' x 5 ' with the constant curvature round metric, and S ^ with the standard 
constant-curvature round metric. 

These wave equations are converted to first-order symmetric-hyperbolic form before solving 
them numerically. The list of dynamical fields u" - {i^, H, O,) is therefore expanded to include 
the first derivatives ofil/.Yl- -dtijj, and <1), - diijj. Constraint damping is used to enforce the 
constraint, 

Ci = - 1>/ = 0, (47) 

using the methods developed in Ref . Il23il with constraint damping parameter - I. 

Exact analytical solutions exist to these wave equations on the three-manifolds used in these 
tests. Therefore the effectiveness and efficiency of the evolution code can be tested in these cases 
by comparing numerical solutions i^/v to this equation with the known analytical solutions ^e- 
The accuracy, and convergence properties, of the code can be measured therefore by monitoring 
the norms of Ai/f - ij/E - 4'n using the solution error measure defined in Eq. ( |33] ). It is also 
useful to monitor the constraint violation errors C,. This is done by constructing the constraint 
error measure: 



This constraint error measure is invariant under changes in the overall scale of the solution, and 
to changes in the coordinates used to represent the solution. 

The tests performed here use the scalar wave evolution system that is implemented as part 
of the SpEC code lEsl E4I1 . This code, developed originally by the Caltech/Cornell numerical 
relativity collaboration, uses pseudo- spectral methods to evaluate spatial derivatives, and the 
method of lines to approximate the hyperbolic system of partial differential equations as sets of 
coupled ordinary differential equations on each collocation point. These tests use an eighth order 
Dormand-Prince Ii25il algorithm to integrate the method of lines ordinary differential equations in 
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time. Each cubic region in these tests is subdivided into one or more computational subregions, 
on which field components are represented using Chebyshev basis functions at the Gauss-Labatto 
collocation points. The structure of these subregions was chosen to achieve fairly uniform spatial 
resolution. The particular choice of subregions is described in the discussion of each particular 
test. 



6.1. Tests of a Multi-Cube Hyperbolic Equation Solver on 

The numerical tests described here use the multi-cube representation of the three-manifold 
with topology given in Appendix A.l The reference metric in this case is the flat Euclidean 



metric, Eq. dA.lb . so the spatial covariant derivatives which appear in the wave Eq. ( |46] | are just 
the Cartesian coordinate partial derivatives. When written in terms of the multi-cube Cartesian 
coordinates on T^, therefore, the wave equation takes the simple form. 



The idea is to solve this equation numerically with initial data 
iJ/(t,x,y,z)\,=Q = cos 
d,if/(t,x,y,z)\,=Q 



2n 

— (kx + iy + mz) 



-w sm 

where k, {, and m are integers, and <jj is given by 

\2 



2n 



{kx + {y + mz) 



The exact solution to this initial value problem is given analytically by 



tl/E(t,x,y,z) - cos 



27T 

(jjt + — {kx + ly + mz) 



(49) 

(50) 
(51) 

(52) 

(53) 



The numerical solution of the wave Eq. ( |49] l for these tests was performed on a set of eight 
computational subregions. These subregions divide the one cubic region needed to represent 
into eight cubes, each half the size of the region in each spatial direction. The internal boundary 
maps between these subregions are just the trivial identity maps. These hyperbolic evolution tests 
were performed using the initial data given in Eqs. dSOl l and ( BTl i with k - £ - m - 2. These tests 
used computational subregions having = 16, 18, 20 and 22 collocation points respectively in 
each spatial direction. The graphs of the solution errors &^ and the constraint violation errors 
&c for these tests are shown in Fig. [T] These graphs demonstrate that the numerical methods 
described here successfully achieve the exponential convergence expected of spectral numerical 
methods. The slow growth in time of the solution error S^, seen in the left side of Fig.|7]is linear 
in time. This type of error is a common feature of the ordinary diff'erential equation integrator 
used for these tests. 



6.2. Tests of a Multi-Cube Hyperbolic Equation Solver on S^ x S^ 

The numerical tests described here use the multi-cube representation of the three-manifold 
with topology 5 ^ x 5 ' given in Appendix A. 2 The reference metric used in this case is the 
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Figure 7: Left: Errors in the numerical solutions hijj tor evolutions witli k = I = m = 2 measured by the 

quantity Right: Constraint errors C, for the evolutions with k = C = m = 2a& measured by the quantity Sq- 



constant-curvature round metric given in terms of angular coordinates [y, 6, ip} in Eq. (IA.8I ). and 
in the multi-cube Cartesian coordinates used in these tests in Eq. iA.9i . This choice of reference 
metric fixes the wave Eq. ( l30l ) to have the form 



■ 5,V + V'V,tA = -5,V + ^ + — + ■, J = 0. (54) 



when expressed in terms of the angular coordinates {x, 0, <p} used on 5 ^ x 5 ' . The idea is to solve 
this equation numerically with initial data: 

il/{t,e,<p,x)r=o = %[e"'^Y(,„(e,'P)], (55) 
drHuO,<p,x)r=o = %[icoe"'^YnM'P)], (56) 

where Y[,„{9, (p) are the standard S~ spherical harmonics, k, {, and m are integers, oj is given by 



Rl rY 



(57) 



and %[Q] denotes the real part of the quantity Q. The exact solution to this initial value problem 
is given analytically by 

Mt^e^f^X) = %[e'^'^""'Ye„,(e,<p)]. (58) 

The numerical solution of Eq. (|54] | is carried out using the Cartesian coordinates of the multi- 
cube description of 5^ x 5 ' described in Appendix A. 2 The spatial covariant derivatives used 



by the SpEC code for this test are evaluated using the Cartesian coordinate representation of the 
round metric given in Eq. ( IA.9I ). The initial data, Eqs. ( l55b and ( l56b . used for these tests are 
evaluated in the multi-cube Cartesian coordinates with the transformations between the angular 
and Cartesian coordinates given in Tables lA~4l and IA3] 

The numerical solution of the scalar wave Eq. ( l54b for these tests was performed on a set 
of twelve computational subregions. These subregions divide the six cubic regions needed to 
represent 5" x 5 cf. Fig. lA.lOl into cubes that are half the size of the region in the z direction. 
The internal boundary maps between these subregions are just the trivial identity maps, while 
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the maps between regions are those given in Table IA.3I These hyperboUc evolution tests were 
performed using the initial data given in Eqs. ( l55b and ( |56] | with k - { - m - 2. These tests 
were performed on computational subregions having = 16, 18, 20 and 22 collocation points 
respectively in each spatial direction. The graphs of the solution errors £^^, and the constraint vio- 
lation errors &c for these tests are shown in Fig. [8] These graphs demonstrate that the numerical 
methods described here successfully achieve the exponential convergence expected of spectral 
numerical methods. The slow growth in time of the solution error seen in left side of Fig.[8]is 
(mostly) linear in time. This growth in the error is a common feature of the ordinary differential 
equation integrator used for these tests. 




t t 

Figure 8: Left: Errors in the numerical solutions hijj for tlie X S ' evolutions with <: = f = m = 2 as measured by the 
quantity £,^. Right: Constraint errors C, for the 5^x5' evolutions with k = £ = m = 2as measured by the quantity &C- 



6.3. Tests of a Multi-Cube Hyperbolic Equation Solver on S ^ 

The numerical tests described here use the multi-cube representation of the three-manifold 
with topology given in Appendix A. 3 The reference metric used in this case is the stan- 



dard constant-curvature round metric for S ^ given in terms of angular coordinates {x, 8, (p] in 
Eq. ( IA.19I ). and in the multi-cube Cartesian coordinates used in these tests in Eq. ( IA.20l i. This 
choice of reference metric fixes the wave Eq. ( |46] | to have the form, 

V7/V7 , d^[sm^xd,ifr] dg [singggiA] . 

Ri^sm-x Ri^ smOsm x Ri^ sm Osm x 

when expressed in terms of the standard angular coordinates {x, 9, ip] used on 5 This equation 
is solved numerically with initial data: 

^{t,9,ip,x),=a = %[YMm(x,e^v)\^ (60) 
d,ilj{t,e,ip,x),=o = %[icjYum(x,S,ip)], (61) 

where YifCm is the S ^ spherical harmonic function defined in |Appendix B[ k, (, and m are integers, 
and to is given by 

(JL> ^ — . (62) 

^3 
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The solution to this initial value problem is given analytically by 



(63) 



The numerical solution of Eq. (|59] | is carried out using the Cartesian coordinates of the multi- 
cube description of 5 ^ described in [Appendix A.3[ The spatial covariant derivatives used by 
the SpEC code for this test are evaluated using the Cartesian coordinate representation of the 
round metric given in Eq. ( IA.20b . The initial data, Eqs. ( l60b and dMT l, used for these tests are 
evaluated in the multi-cube Cartesian coordinates with the transformations between the angular 
and Cartesian coordinates given in Table lA.8] and lA.9l 

The numerical solution of the scalar wave Eq. ( |59] ) for these tests was performed on a set 
of eight computational subregions. These subregions are identical to the eight cubic regions 
needed to represent S^, cf. Fig. lA.l II The maps between regions are those given in Table lA.Tl 
The hyperbolic evolution test was performed using the initial data given in Eqs. (l60t and dHTT l 
with k - ( - m - 2. These tests were performed on computational subregions having = 
16, 18, 20 and 22 collocation points respectively in each spatial direction. The graphs of the 
solution errors &^ and the constraint violation errors &c for these tests are shown in Fig. |9l 
These graphs demonstrate that the numerical methods described here successfully achieve the 
exponential convergence expected of spectral numerical methods. The slow growth in time of 
the solution error seen in the left side of Fig.|9]is (mostly) linear in time. This growth in the 
error is a common feature of the ordinary differential equation integrator used for these tests. 
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Figure 9: Left: Errors in the numerical solutions Ai/< for the 5 evolutions with k = C 



quantity £^ . Right: Constraint errors C, for the 5 ' evolutions with k ■ 



2 as measured by the 



: 2 as measured by the quantity &c- 
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Appendix A. Examples of Multi-Cube Representations of Three -Manifolds 

This appendix describes the construction of multi-cube representations of manifolds using 
the methods developed in Sees. |2] and [3] Each multi-cube representation consists of a set of 
non-overlapping cubes Sa that cover the manifold, a set of maps ^F^^ that identify the faces 
of neighboring cubes, and finally a smooth positive definite reference metric gij used to define 
the differential structure on the manifold. The construction of these multi-cube structures is 
described here for three common three-manifolds: the three-torus with a flat reference metric, 
the spherical-torus 5' x 5 ' with a constant-curvature round-sphere metric, and the three-sphere 
S ^ with the standard constant-curvature round-sphere metric. These examples are used in Secs.|5] 
and|6]to illustrate the solution of partial differential equations on multi-cube manifolds using the 
methods developed in Sec. |4] 

Appendix A.l. Multi-Cube Representation of 

The simplest example of a multi-cube manifold is the three-torus, . Only a single cube B\ 
is needed to cover this manifold, and it is most convenient to locate this cube at the origin in 
so ci = (0, 0, 0). Opposite faces of this cube are identified without rotation or reflection to obtain 
the topology: « d-^S^ d+ySi ^ d-ySi, and d+^Si ^ d-,!Bi. The maps, >fJ^-;, 

and ^Pj^!, needed to effect these identifications are defined by Eq. ([1]) with the rotation matrices, 

Cg^, being just the identity matrices: C[;^'J = C[^J, = C[^: - I. The three-torus admits a 
smooth flat metric, so a convenient choice of reference metric for this manifold is: 

ds~ - gijdx'dx-' - dx^ + dy^ + dz^, (A. 1) 

where x, y and z are the multi-cube Cartesian coordinates that label points in Si . 

Appendix A.l. Multi-Cube Representation ofS~ X S ' 

The manifold 5 ~ x 5 ' can be covered by a set of six cubic regions Sa with a - {!,..., 6). A 
convenient way to arrange these cubes in R^ is illustrated in Fig. lA.lOl The values of the cube- 
center location vectors ca for this configuration is summarized in Table I A. 2l The inner faces of 
the touching cubes in Fig. lA.lOl are connected by identity maps, while the outer faces are identi- 
fied using the maps described by Eq. ([T]) with the rotation matrices C^^ given in Table lA.3] This 
representation of 5^ x 5 ' is constructed by taking the Cartesian product of S ' (the periodically 
identified z-axis in this representation) with the commonly used "cubed-sphere" representation 

of^^ciaa. 



Table A. 2: Cube-Center Locations for 5" X S ' 



?1 


= (0, -L,0) 




= (0,L,0) 


C5 


= (i,0,0) 


C2 


= (0,0,0) 


C4 


= (0,2L,0) 


C6 


= (-L,0,0) 



It is useful to discuss the method used to construct the "cubed-sphere" representation of 5 ^ in 
some detail here, since this method is used in Appendix A. 3 as the model for constructing a new 
representation of 5 . Let {x,y,z] denote Cartesian coordinates in an R^ , and let -H -H = r^ 
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Figure A. 10: The three-manifold 5^ X S ' is represented using the six cubic regions illustrated here. The faces of these 
cubes are identified using the maps described in Table Ia31 This representation of X S ' is based on the commonly 
used "cubed-sphere" representation of 5^. 



Table A. 3: Cube Face Identifications. (9„Sa (i/sSs. and rotation matrices. CgS, for the interface maps in 5^ X S ' . 









r^Aa 
^8/3 


'-'At? 




«-> 




r^Aa 


^Aa 




«-> 




I 


I 




<-> 




I 


I 




«-> 




I 


I 




<-> 




R+ 


R 




«-> 




R , 


R+ 




<-> 




I 


I 




<-> 


d-yS3 


I 


I 




<-> 




I 


I 




<-> 




I 


I 




<-> 




I 


I 




<-> 


5-,S4 


I 


I 




<-> 


d^ySi 


R , 


R« 




<-> 




R+ 


R c 




«-> 




I 


I 




«-> 




Rl 


Rl 


5-.S4 


«-> 


d-,S6 


Rl 






«-> 




I 


I 


5+cS6 


<-» 




I 


I 



denote a two-sphere 5 ~ of radius r. It is useful for some purposes to identify points on this S ^ 
using standard angular coordinates and tf. 

X - r sin 6 cos (fi, (A. 2) 

y - rsinOsiiKf, (A. 3) 

z - rcosO. (A.4) 

Now consider a cube S centered at the origin, of size L - 2r/ V3 (which just fits inside the 
sphere), whose orientation is aligned with the {x,y,z} axes. Let daS represent the six faces of 
this cube, with a = +x, etc., labeling the various faces. The images of these six faces can be 
arranged in a plane, like the a = +z faces of the cubes shown in Fig. lA.lOl The goal here is to 
construct a representation of 5 " x 5 ' , so it will also be useful to make a correspondence between 
these cube faces dcrS with the cubes shown in Fig. lA.lOl Table lA.4l gives the relationship between 
the cube-face identifiers a - +x, etc. and the cubic region labels a=i,2,...,6 shown in Fig. lA.lOl 
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Points on each of the cube-faces, dsS, can be identified by their local Cartesian coordinates. 
For example, points on the a = +z face, i.e. the a = 2 face in Fig. I A. 101 can be identified by 
the coordinates {x,y}. It is also useful to introduce scaled local Cartesian coordinates, {Xa, Ya} 
to represent the points on these faces. For the a = +| face for example, it is useful to set 
{X2,Y2} = {x/z,y/z}- Each coordinate has been divided by z, which is constant on this face, 
to ensure that the scaled coordinates {X2, Fa) are confined to the ranges, -1 < X2 < 1 and 
-I < Y2 < I. Similar definitions are made on the other faces, cf. Table lA~4l that ensure 
the Xa and Ya are all oriented the same way as in Fig. lA.lOl and all satisfy -I < Xa < I and 
-I < Ya < I - Using Eqs. ( IA.2l )-( IA.4b . this construction provides a natural identification between 
points on the original sphere, labeled by their angular coordinates {0, (p), and the Cartesian cube- 
face coordinates {Xa, Ya} via the equations summarized in Tables lA4l and lA3] 



Table A. 4: Cubed-Sphere Representation of S " : Angular to Cartesian Coordinate Map. 



A 


a 


Xa 


Ya 


1 


-y 


X _ 

V 


— COt(f 


y 




- cot CSC (f 


2 


+z 


A _ 


tanflcosi/j 


y 




tan sin (f 


3 


+y 


X _ 

y ~ 


COtl^ 


y 




- cot CSC (f 


4 


-z 


_x _ 


- tan cos (f 


y 




tan sin (f 


5 


+x 


X 


-cot0sec(/) 


V 
X 




tan(/? 


6 


—X 


X 


- col stop 


y 

X 




-tan(/3 



The [Xa, Ya] defined in this way are local Cartesian coordinates. These could be converted 
to global coordinates by adding in the appropriate offset for each face: jc^ = -1- ^LXa and 

= c\ + ^LYa- Alternatively, the angles tan"' Xa and tan"' Ya could be used as local "Cartesian" 
coordinates on these cube faces. These angle-based Cartesian coordinates have the advantage 
of giving a more uniform mapping of the Euclidean plane onto the image of the cube face on 
the sphere, so they are the preferred choice for numerical work. Global Cartesian coordinates 
constructed from these angle-based coordinates are defined by 

2L 

x\^c\ + — tan ' Xa, (A.5) 
n 

4 = < + — tan"' Ya, (A.6) 
n 

where Xa and Ya are functions of the standard angular coordinates and ip by the expressions 
given in Table lA.4l 

For representations of 5 ^ x 5 ' , an appropriate coordinate is also needed for the periodically 
identified z direction in Fig. lA.lOl Introduce an wgle whose range is —n < x < t^, that labels 
the points in the S ' subspace. Then define the global Cartesian coordinate associated with this 
direction as 

4 = 4 + ^X, (A.7) 
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Table A.5 : Cartesian to Angular Coordinate Map for the Cubed-Sphere Representation of S ^ . The range of the local 
Cailesian coordinate is - 1 < Xa < 1 , and the range of (? is < 6 < in these expressions. The ranges of <f for 
different values of Ya are specified in the table. 



-range 



-1 < Fi < 1 

1 > ^2 > 
-1 < ^2 < 

-1 < ^3 < 1 
1 > ^4 > 
-1 < ^4 < 

-1 < Fs < 
1 > ^5 > 

-1 < ^6 < 
1 > ^6 > 



COS (fi 



Xi/^Xl + Y^^ 
X^/^l+Xj 
XJ ^Xj + Yl 
X,l ^Xl + Yl 



II Jl + Yl 



i^-range 



Ik 



5n 



ji>(f>Q 
2j: > ip > n 

2n > if > n 
n>(p>0 

2n>ip>^-f 
f >^>0 



5;r 



> if > n 



COS 



Yd^l+X\ + Yl 
\I^\+Xl + Yl 
l/^l+Xj + Yl 
-Yi/^l+Xj + Yj 
-l/^l+Xl + Yl 
-l/^l+Xl + Yl 



-Xs/^l+Xl + Yl 
-X,l ^l+Xl + Yl 
Xd^\+Xl + Yl 
Xe/yjl+Xj + Yl 



The standard constant-curvature "round" metric on 5^ x 5 ' is smooth, and it is therefore an 
acceptable choice for the reference metric to define the diff'erential structure on this manifold. 
The simplest representation of this round metric uses the angular coordinates 0, (p, and;;^: 



ds- = Rl{d0~ + sin' 0d(p-) + R^dx', 



(A.8) 



where R2 and Ri are constants that specify the radii of the and S ' parts of the geometry re- 
spectively. Using the transformations given in Eqs. (IA.5l )-( lA77l i and Table lA.4l a straightforward 
(but lengthy) calculation gives the global multi-cube Cartesian-coordinate representation of this 
metric on 5 ^ x 5 ' : 



ds-" 



(nR2\- il+Xl){l + Yl) 
(1+Xl + Yj)2 



lnR2 \- 
l 277/ 



[(1 + Xl){dxlf - IX^YAdx^d^^ + (1 + Yl){d^^f] 



IttR, 



idr.f 



(A.9) 



The Xa and Ya that appear in this expression are thought of as the functions of the Cartesian 
coordinates obtained by inverting the expressions given in Eqs. ( IA.5I ) and ( IA.6I 1: 



Xa - tan 
Ya - tan 



2L 



2L 



(A. 10) 
(A.ll) 
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The functions Xa and Ya depend on the location of a particular coordinate region through the pa- 
rameters c;^ and c\. However, beyond this dependence the multi-cube coordinate representation 
of the 5^x5' round metric given in Eq. iA.9i is the same in each of the six coordinate regions 

These multi-cube Cartesian coordinates {xA,yA,ZA} turn out to be harmonic with respect to 
the round metric on 5^ x 5 i.e., each coordinate is a solution (locally within each cubic-region, 
not globally across the interface boundaries) to the covariant Laplace equation, - V^^, xa - 

Vai yA - ^'a ^Ai za , where V/i,- is the covariant derivative associated with the 5^x5' metric in 
region a. These conditions are equivalent to = Bai ( y/gAgA) where gA - detgAij and g'j is the 
inverse of the metric gAij expressed in terms of the multi-cube Cartesian coordinates in region a. 

Appendix A.3. Multi-Cube Representation of 

The locations of the eight cubic regions used to construct this representation of are illus- 
trated in Fig. lA.llI The values of the cube-center location vectors ca for this configuration is 
summarized in Table IaTSI The inner faces of the touching cubes in Fig. lA. Ill are assumed to 
be connected by identity maps. The outer faces of these eight cubic regions are identified using 
the maps described in Table IA.7I This "cubed-sphere" representation of 5 ^ is a natural three- 
dimensional generalization of the two-dimensional cubed-sphere representation of described 
in Appendix A. 2 It is constructed by inserting a four-dimensional cube into a three-dimensional 
sphere 5"* inR*, and then identifying points on the faces of the four-cube with the points on the 
three-sphere that are connected by rays extending outward from the origin. 




Figure A. 11: The three-manifold 5^ can be represented using the eight cubic regions illustrated here. Cubic region S2, 
centered at the origin (?2 = (0, 0, 0) is hidden between S^ and Sg in this figure. The outer faces of these cubes are 
identified using the maps described in Table lA.Tl 



Table A. 6: Cube-Center Locations for 





= (0,-L,0) 


C3 


= (0,L,0) 


C5 = (L, 0, 0) 


c^ = (0,0, L) 


Q 


= (0,0,0) 


C4 


= (0,2L,0) 


C6 = (-L,0,0) 


= (0,0,-L) 
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Table A. 7: Cubic Region Face Identifications, daSA <-» O/iSb , and rotation matrices, Cgg, fortlie interface maps in 5^. 
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^8/3 




5„2 
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B 


C^Aa 
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R , 




5i 


<-> 
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R . 
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d-A 
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I 
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d^ri 


S3 


<-> 


d-y'l 
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<-> 
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R+ 


d-J 


53 


<-> 


d^yi 
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R+ 


R - 




<-> 
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d-J 


53 
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d+yl 
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R V 


R+.r 




«-> 




R;. 


R;. 


d-J 


34 


«-» 




R;. 


Rl 




<-> 




R;. 


Ri. 




<-> 


d-zi 
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R;. 


Rl. 




<-> 




Ry 


R+v 


d-zi 


55 


<-> 


d^,i 
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R+v 


Ry 




<-> 




R+v 


R y 


d-zi 


56 


<-> 




R y 


R+v 



It is appropriate to discuss this "cubed-sphere" representation of 5 in some detail, since it 
does not appear to have been used or described in the literature before. Let {x,y,z,w} denote 
Cartesian coordinates in R'^, and let + z' + w' = denote a three-sphere, 5'"', of radius r. 

It is often useful to identify points in 5^ using the angular coordinates 9 and (p: 

X - r sin^ sin 9 cos ip, (A. 12) 

y - r smx sin 9 sin ip, (A. 13) 

z = r sin;^' cos 0, (A. 14) 

M> - rcos;if. (A. 15) 

Now consider a four-cube centered at the origin, of size L - r (which just fits inside the three- 
sphere), whose orientation is aligned with the {x,y,z,w] axes. Let daS denote the eight faces 
of this four-cube (each of which is a three-cube) labeled by the index a - +x, etc. Arrange 
the images of these eight three-cubes in R? at the locations given in Table IA.6I as shown in 
Fig. lA. Ill Table IATS] gives the relationship between the four-cube face identifiers a - +x, etc. 
and the three-cube region identifiers a=i,2,...,8 shown in Fig. lA. Ill 

Points on each of the four-cube faces, daS, can be identified by their local Cartesian coor- 
dinates. For example, points on the a - +w face, i.e. the a - 2 region in Fig. lA.l II can be 
identified by the coordinates {x,y,z]. It is convenient to introduce scaled local Cartesian coordi- 
nates, {Xa, Ya,Za} to represent the points on these faces. For the a = +w face for example, set 
{X2, Y2,Z2] = {x/w,y/w,z/w}- Each coordinate has been divided by w, which is constant on this 
face, to ensure that the scaled coordinates {X2, Y2,Z2} are confined to the ranges, -1 < X2 < 1, 
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< Y2 < 1, and -1 < Z2 < 1. Similar definitions are made on the other faces, cf. Table I A. 8l 
that ensure the Xa, Ya, and Za are all oriented the same way as in Fig. lA.llI and all satisfy 
-1 < Xa < 1,-1 < Ya < 1, and -1 < Za < 1. Using Eqs. (IA.12l i-( IA.15l l, this construction pro- 
vides a natural identification between points on the original three-sphere, labeled by their angular 
coordinates {x, 0, ip}, and the local Cartesian coordinates {Xa, Ya, Za} on each four-cube face via 
the equations summarized in Tables lAT8] and lA!9] 



Table A. 8: Cubed-Sphere Representation of 5^. 
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X 


- - cot CSC sec ip 




- 1 = - tan If 
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= - cot 9 sec (fi 
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1 - tan cos 1^ 




\ - tan 9 sin ip 
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- - cotx sec 9 
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-4 = - tan cos (p 




- ^ = - tan 9 sin ip 
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- - cot;^' sec 9 



The {Xa, Ya,Za} defined using this cubed-sphere construction are local Cartesian coordinates 
on each of the faces of the four-cube. They could be converted to global coordinates by adding 
the appropriate offset for each cube: x-^ - c'^ + \LXa, x^^ - + \LYa, and - c~^ + \LZa- 
Alternatively, the angles tan ' Xa, tan ' Ya, and tan ' Za also provide local Cartesian-like coor- 
dinates for these cubes. These angle-based Cartesian coordinates give a more uniform mapping 
of Euclidean space onto the image of the four-cube face on the three-sphere. So as in the two- 
dimensional cubed-sphere case, these angle-based Cartesian coordinates are the preferred choice 
for numerical work on the multi-cube representation of 5^. Global multi-cube Cartesian coordi- 
nates constructed from these angle-based coordinates are defined by 

4 = + — tan"' Xa, (A. 16) 

TT 

= + — tan"' i'^' (A- 17) 

n 

yr. ^c'.+— tan"^ Za, (A. 18) 

n 

where Xa, Ya, and Za are functions of the hyper-spherical angular coordinates 9 and ip given 
by the expressions in Tables IAT8] and |A!9l 

The standard constant-curvature "round" metric on 5 ^ is smooth, and it is therefore an ac- 
ceptable choice for the reference metric to define the differential structure on this manifold. The 
simplest representation of this round metric uses the angular coordinates x, 6, and p: 

ds^ = Rl{dx^ + sm^Xd9'^ + ^'^^^X^'^^'Sdip^)' (A.19) 
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Table A. 9: Cartesian to Angular Coordinate Map for the Cubed-Sphere Representation of S'. The range of the local 
Cartesian coordinate A"^ is -1 < < 1, the range of is -1 < Za < 1, the range of the angular coordinate 6 is 
< < n, and the range of ^ is < ^ < in these expressions. The ranges of (f corresponding to different ranges of Ya 

are specified in the table. The quantities Wa = -^i + + are used to simplify the expressions for cos^. 



Fa -range 



-1 < Fi < 1 

1 > ^2 > 
-1 < ^2 < 
-1 < 73 < 1 

1 > 74 > 

-1 < ^4 < 

-1 < Fs < 
1 > > 
-1 < ^6 < 

1 > ^6 > 

1 > Fv > 
-1 < Fv < 

1 > Fs > 
-1 < Fs < 



cos (f 



Xl + F| 



X.J^Xl + Yl 
X^l^l+X] 
X^l ^Xl + Yl 
X,l ^Xj + Yj 

XtI ^X^ + f2 
Xt/ ^X'- + Y] 
Xdyjxl + Yl 
Xf.1 ^Xl + Yl 



i/j-range 



TT > 1^ > 

2n > ip>Ti 

2j: > ip > n 
n>ip>Q 

2n>ip>^-f 
f >^>0 

f ^^>^ 

^>^> f 

TT > 1^ > 

2n > if > n 
n>ip>Q 
2j: > ip > n 



cos 9 



Zi/^1 +x\+z\ 



xl + Yl + Zl 



ZoJ^Xl + Yl+Zl 

Z3/^i +xl + zl 



Xl + Yl + Zl 



Z4/^ 

Z4/ ^Xl + Yl + Zl 
Zs/^l+Yl+Zl 
Zs/^l + Yl+Zl 
ZJ^l + Yl+Zl 
ZJ^l+Yl+Zl 
l/^l+Xl + Yl 
l/^l+Xl + Yl 
-l/ ^l+Xl + Yl 
-l/^l+Xl + Yl 



cosx 



I/W2 
I/W2 
-Yi/Wi 

-im 
-im 
-X5/W5 
-X5/W5 

Xem 
Xem 
-Zt/Wi 
-Zt/Wi 

Zg/Wg 



where is a constant that specifies the radius of the S^. Using the transformations given in 
Eqs. (IA.16b -( lA.18b and in Tables IATS] and |A!9l a straightforward (but lengthy) calculation gives 
the global multi-cube Cartesian-coordinate representation of this metric on 5^: 



IKR,\ HI+Xl)(l + Yl){l+Zl) 
\2l) (l+xl + Yl+Zl)^ 



{l+Xl)(l+Yl+Zl) 



(1+F2)(1+Z2) 



idxiy- - 



(l + Yl){l+Xl+Zl ), 
(1+X2)(1+Z2) 



(1+Zl)(l+Xl + Yl) 

-(dr.) - 



2XaYa 

i+zl 



l+Yl ^ 

2YaZa V , - 
1 + x^ "^a"^a 



dxl dx\ 



(A.20) 



(1+Xl){l + Yl) 

The Xa, Ya, and Za that appear in Eq. ( IA.20b are thought of as the functions of the global multi- 
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cube Cartesian coordinates obtained by inverting the expressions given in Eqs. ( IA.16b - (IA.18b : 



Ya 
Za 



tan 



tan 



tan 



2L 

/ V V \ 



IL 

^(4 - 4) 



IL 



(A.21) 
(A.22) 
(A.23) 



The functions X^, Ya and Za depend on the location of a particular coordinate region through 
the parameters c;^, and c^. However, beyond this dependence the multi-cube coordinate rep- 
resentation of the S ^ round-sphere metric given in Eq. ( IA.20b is the same in each of the eight 
coordinate regions 'Ba- 

These multi-cube Cartesian coordinates \xA,yA,ZA\ turn out to be harmonic with respect to 
the round metric on 5^, i.e., each coordinate is a solution (locally within each cubic-region, not 
globally across the interface boundaries) to the covariant Laplace equation, = '^'a^ m xa - 
^A^AiyA = ^'a^AiZa, where V^; is the covariant derivative associated with the metric in 
region a. These conditions are equivalent to = 5a, ( -\/gAg'J) where gA - detgAij and g'j is the 
inverse of the metric gAij expressed in terms of the multi-cube Cartesian coordinates in region a. 



Appendix B. Spherical Harmonics on 

This appendix derives expressions for the eigenfunctions of the Laplace operator on the three- 
sphere S^. These eigenfunctions are referred to here as three-sphere harmonics. These functions 
are defined as solutions of the equation 

Y'^iY^-AY, (B.l) 

where V, is the covariant derivative operator on 5^, and A is an eigenvalue. These functions have 
been studied previously by a number of authors 1 2^ 27 , 28 , 2^ . Here a slightly different repre- 
sentation is introduced that allows these harmonics (of arbitrary order) to be evaluated accurately 
in a straightforward way. Using the angular coordinate representation of the round metric on S ^ 
from Eq. ( IA.19I ), it is straightforward to write the co-variant Laplace operator explicitly as 

Ri^sirr X 7?^ sin 6? sin" Ri^sirr 6sm x 



The eigenvalue problem, Eq. ( IB. lb . can be solved then by separation of variables. The non- 
singular solutions to this equation have the form: 

Nil- f+ 1 

YkCm(x, 9, v) = QJ (cos;^)r/(cos 0)e""^ (B.3) 

where Py and ^ are the associated Legendre functions of the first and second kind respectively. 
The eigenvalue associated with this Kj^f,,, is 

k(k + 2) 
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These functions are non-singular on 5 ^ only for integers k, ( and m satisfying 



k>Q, 
k>{>0, 
(>m> ~l. 



(B.5) 
(B.6) 
(B.7) 



The half-integer associated Legendre functions Q^^\ W with x - o.o'&x are non-singular for 
- 1 < X < 1 , and can be evaluated re-cursively. For fixed I, the functions with A: < ^ can be shown 
to vanish. 



(B.8) 



using §3.4 Eq. (13) in Ref. ll30ll . For k - I similar argument using §3.6.1 Eq. (14) in Ref. |13; 
gives 



d;}M) =(-i)^^'2^^!Vf(i--^T 



(B.9) 



The functions with k> I can be determined from these using the recursion relation, 

{k-i+ 2)eJ^J(x) = i(k +2)x e[^|(x) -{k+e+ 2)d^}M), 



from §3.8 Eq. (12) in Ref. tM- Evaluating Eq. dRTOt for ;t = ^ - 1 gives 

d^^{x)^2(l+\)xd^}^{x). 



(B.IO) 



(B.ll) 



using Eq. ( IB. 8b . The Q^^\ (x) with k > { + 2 can then be generated recursively using Eq. ( IB. 10b . 

This recursion relation is known to be a stable and accurate way to generate the Legendre func- 
tions of the first kind, P"\x), cf. Ref. 113 ill . Our numerical tests indicate that it is also an accurate 

way to generate the half-integer Legendre functions of the second kind, Q ] (x). 

2 

The orthogonality properties of the Ykimix, 0, ip) are determined by the orthogonality proper- 
ties of 1 (cos;if), P^,(cos 9) and e'""^. The needed condition for 2^^^ can be obtained from the 
associated Legendre differential equation. 



dx 

from which it follows that 



(1 )—r- 

dx 



v(v +\)- 



1 



(B.12) 



d 
dx 



dx 



{V -v){v + v' + \)^.Q^',. 



(B.13) 



The half-integer associated Legendre functions are well behaved in the interval -1 < x < 1, 
therefore integrating Eq. ( IB.13I ) over this interval gives 



= (v' - v)(v + v' -I- 1) 



x)dx. 



(B.14) 
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l+- 

It follows that the , (x) with k>0 and t >Q satisfy the orthogonality condition: 



where M^t is the numerical constant. 



Mi 



nHk + (+\)\ 



'''^ 4(k+i)(k-£y.' 

The analogous orthogonality relations for P^/cos 6) and e'""^ are well known: 



where 



1.71 6 m'm 



Jo 



■"'(y)Pf(y)dy, 



(( + m)\ 



From these conditions then, it follows that by choosing the normalization constants 

1 

Nktm = — ^ , 

the Yu„, satisfy the following orthogonality conditions on 5^, 

2 . 



(B.15) 
(B.16) 

(B.17) 
(B.18) 

(B.19) 
(B.20) 



Yu.v„,Yli„^^d'^x = r\\ dx\ de\ d^Wx^vciQYvfniYl 

Jo Jo Jo 



27r Jo 



(B.21) 
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